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We present calculations of the ground state and excitations of an anisotropic dipolar Bose gas 
in two dimensions, realized by a non-perpendicular polarization with respect to the system plane. 
For sufficiently high density an increase of the polarization angle leads to a density instability of 
the gas phase in the direction where the anisotropic interaction is strongest. Using a dynamic 
many-body theory, we calculate the dynamic structure function in the gas phase which shows the 
anisotropic dispersion of the excitations. We find that the energy of roton excitations in the strongly 
interacting direction decreases with increasing polarization angle and almost vanishes close to the 
instability. Exact path integral ground state Monte Carlo simulations show that this instability 
is indeed a quantum phase transition to a stripe phase, characterized by long-range order in the 
strongly interacting direction. 



Strongly correlated dipolar Bose gases in two dimen- 
sions (2D) polarized along the direction normal to the 
system plane have been extensively investigated in re- 
cent years QUI]. The ratio between the dipolar length 
ro = mCddl '(^Trh 2 ) and the average interparticle distance 
provides a measure of the strength of the interaction. Cdd 
is the coupling constant proportional to the square of the 
(magnetic /x or electric d) dipole moment, resulting in a 
dipolar length that can range from a few A for mag- 
netic dipolar systems like 52 Cr (/1 = Q^b, with /i^ the 
Bohr magneton) , to thousands of A for heteronuclear po- 
lar molecules like KRb, LiCs [5], or RbCs [6 . However, 
chemical reactions and three-body losses impose limita- 
tions on what can be measured in experiments with po- 
lar molecules. Therefore, recent efforts focus also on ex- 
otic lanthanide magnetic systems like 164 Dy or 168 Er, [7 
where the combined effect of a large magnetic moment 
(/i = lOfiB for 164 Dy and /1 = 7/jb for 168 Er) and a large 
mass, lead to dipolar length scales that, although still 
significantly lower than the corresponding value for po- 
lar molecules, is several times larger than that of 52 Cr. 
Er2 with fi = and twice the mass of Er would reach 

even higher values of [8]. 

A 2D dipolar Bose gas polarized along the normal 
direction to the confining plane develops a roton exci- 
tation at high density due to the strong repulsion be- 
tween dipoles at short distances [3]. Other works have 
revealed competing effects in a quasi-2D geometry due 
to the head-to-tail attraction of the dipole-dipole inter- 
action when the third spatial dimension is added, to the 
point that the system becomes unstable against density 
fluctuation below a critical trapping frequency in that 
direction [5HTT]. This leads to the question of whether 
a similar situation can hold in a purely 2D geometry 
when a head-to-tail component to the dipole-dipole in- 
teraction is added by tilting the polarization with respect 
to the direction normal to the system plane. The inter- 
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, with particles moving in the x, y-plane and 
a polarization field in the x, z-plane, tilted by an angle 
a with respect to the z-axis. The interaction is weak- 
ened in the x-direction as a is increased, while it does 
not change in the y-direction. Notice that, in the case of 
bosonic particles, only polarization angles where V(r) is 
non-negative, i.e. a < a c = arcsin(l/v / 3) = 0.61548 . . ., 
are meaningful, if there is no additional short-range re- 
pulsion to prevent two dipoles from collapsing to a point. 

The effect of a tilted polarization on the superfluid 
response of a quasi-2D dipolar Bose gas has been inves- 
tigated by mean field theory [IT]. The appearance of 
a stripe phase has been predicted in 2D dipolar Fermi 
systems by approximate methods [T^HH], observed as 
a spontaneous symmetry breaking even in the isotropic 
case (a = 0) for high interaction strength. However, re- 
cent fixed node diffusion Monte Carlo simulations find no 
evidence of that in the isotropic case before the system 
crystallizes [15]. In previous work [16] we investigated 
the low density regime of the 2D dipolar Bose gas of 
particles interacting by the potential U(r), analyzing the 
universal energy scaling properties of the anisotropic gas 
and other ground state properties. Up to the maximally 
allowed polarization angle a c , the low-density system al- 
ways remains in a gaseous form and no trace of a stripe 
phase is found. In this work we focus on the high den- 
sity regimes of this system, studying the effect of the 
anisotropy on the dispersion relation, especially the ro- 
ton, and show how a stripe phase forms at large densities 
and polarization angles. Throughout the paper, lengths 
and energies are given in units of r and E = 7i 2 /(mr 2 ), 
respectively. 

Before showing exact ground state results obtained 
by Monte Carlo simulations, we present a qualitative 
stability analysis of the ground state for a wide range 
of densities n and polarization angles a. We use the 
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hyper-netted chain Euler-Lagrange (HNC-EL) [T5] tech- 
nique, which is based on a Jastrow-Feenberg ansatz for 
the bosonic many-body wave function $(ri, . . . , rjv) = 
expJ2i<j u 2( r i — r j)- We determine u 2 (r) variationally 
by imposing the extremal condition . = for the 



ground state energy E g and solving for the pair distri- 
bution function g(r) within the approximate HNC-EL/0 
framework. It was shown in Ref. [19] that when the low- 

est eigenvalue Aq of the hessian K(r, r ) 



is non-positive, the system becomes unstable against 
infinitesimal fluctuations of g(r), where the associated 
eigenfunction /o (r) is the fluctuation driving the instabil- 
ity. We have obtained the lowest eigenvalue and eigenvec- 
tor of K(r : r') by imaginary-time propagation. Figure [l] 
shows Ao(a) as function of a for a wide range of densities, 
where each curve is normalized by the respective isotropic 
limit, Ao(0). For all n, Ao(a) decreases with increasing a, 
but there is an important distinction between its behav- 
ior at low and high densities: for n < 128, Xo(a) remains 
finite up to a c , while for n > 128 Ao(a) falls to zero 
already before reaching a c . Hence the high-density gas 
state is unstable above a critical angle ao which is smaller 
than a c ; this is also seen by the fact that the HNC-EL 
equations do not converge in the range ao < a < a c . The 
inset of Fig. [I] shows fo (r) for n = 256 at the largest an- 
gle where we found solutions to the HNC-EL equations, 
a = 0.58, where A o (a o )/A o (0) almost vanishes. Away 
from the correlation hole at r = 0, /o( r ) is essentially a 
plane wave in the more repulsive y-direction. indicating 
a tendency of <?(r) towards long-range order in the y- 
direction. We stress that the HNC-EL/0 results become 
less accurate with larger n, hence the stability limits are 
only approximate. Exact simulation results are presented 
below. 

The anisotropic nature of the interaction which desta- 
bilizes the system beyond ao also influences the spectrum 
of elementary excitations which we investigate by calcu- 
lating the dynamic structure function S(k, E). 5(k, E) is 
proportional to the probability that a perturbation trans- 
fers momentum k and energy E to the system. Thus, for 
a given k, S(k, E) has a marked peak if E coincides with 
the energy of an excitation of the system. We obtain 
5(k, E) using the dynamic many-body theory |20j . where 
the equations of motion for time-dependent fluctuations 
of up to pair-correlations in the many-body wave function 
are solved numerically. If the convolution approxima- 



21] is used 



tion for the three-body distribution function 

S(k,E) is obtained as S^k, E) = — g _ S ( k E ^ 

where Y,(k,E) is the complex, energy-dependent self- 
energy of Eq. (2.46) in Ref. [2D]. We note that the only 
input required to calculate £(k, E) is the static structure 
factor S(k) of the ground state. 

In order to get exact results for 5(k), we have carried 
out stochastic simulations using the path integral ground 




FIG. 1. The lowest eigenvalue Ao(a) of the hessian of the 
ground state energy E g is shown in HNC-EL/0 approxima- 
tion for densities n — 8; 32; 128; 256. Ao(a) is normalized by 
the respective eigenvalue at a — 0. The inset shows the eigen- 
function /o(r) for n = 256 and a = 0.58. 



state (PIGS) Monte Carlo technique of Ref. [22] which 
starts from a variational wave function </>o and projects 
out components orthogonal to the true ground state by 
propagation in imaginary time. In this sense, the result 
of the simulation becomes stochastically exact provided 
the approximation employed for the Green's function is 
accurate and the propagation time is long enough [23] . In 
the present case we have used as (f>o a Jastrow-Feenberg 
ansatz Yii<j f{ r ij) built from the two-body correlation 
factor /(r) = K (2/^/r), corresponding to the zero- 
energy solution of the two-body problem of the isotropic 
1 /r 3 interaction. Despite the isotropy of tfio , anisotropic 
contributions are taken into account by a fourth-order 
propagator [24] , which contains the anisotropic potential 
V(r) and its gradient. The results presented in this work 
have been obtained for N = 512 particles in a simula- 
tion box with periodic boundary conditions. Addition- 
ally, simulations with smaller N have been carried out 
in order to see the A^-dependence of the highest peaks in 
S(k) when n and a increase. 

In Ref. [3] we studied the density dependence of 5(k, E) 
of the 2D dipolar quantum gas in the isotropic limit (a = 
0). There we found a spectrum with a pronounced roton 
for large density, due to the strong correlations induced 
by the 1/r 3 repulsion. Here, we are interested in the 
dependence of 5(k, E) on the polarization angle a and, 
for a given a > 0, its dependence on the direction of 
k. In Fig. [2] we show S*(k, E) for n = 128 and a = 
0.20; 0.50; 0.58 in order to illustrate the evolution from 
an isotropic to an anisotropic excitation spectrum and 
the approach to the stability limit. The wave vector k 
is pointing in the y and x-direction (i.e. the direction of 
strongest and weakest interaction) in the left and right 
panels. We broaden 5(k, E) by adding a small imaginary 
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FIG. 2. S(k, £) for k = (0, ft) (left panels) and (ft, 0) (right 
panels) for polarization angles a = 0.20; 0.50; 0.58 at density 
n = 128. The spectrum in Bijl-Feynman approximation is 
shown as a solid line, and the dotted line denotes the damping 
limit E e (k). 



part rj = 0.2 to the energy in the calculation of S(k, E), 
since otherwise undamped modes would not be visible in 
Fig. [2j Also shown is the Bijl-Feynman approximation 
of the spectrum, obtained by setting S(k, E) = (solid 
line) . 

For a = 0.20 the dispersion is almost independent 
on the direction of k, with only a slight slope of the 
Pitaevskii plateau [25], which for isotropic systems de- 
notes the sudden onset of damping at twice the roton 
energy due to decay into two rotons. As a is increased, 
S{k, E) becomes very different in the y- and x-direction 
and features a highly anisotropic dispersion relation for 
a = 0.58. The wave number of the roton depends on 
the direction of k, but most strikingly its energy decays 
almost to zero in the y-direction for a = 0.58, indicating 
that the system is close to the limit where the homo- 
geneous gas phase in unstable against infinitesimal den- 



sity fluctuations. Since the restriction to pair correlation 
fluctuations used here typically gives an upper bound to 
the excitation energy [26], the exact roton energy in y- 
direction is expected to be even smaller. Furthermore, 
at twice the wave number of the roton, S(k, E) has an- 
other roton-like peak for a = 0.58, following a quadratic 
dispersion, albeit broadened and with smaller spectral 
weight. In the y-direction, the dispersion relation thus 
resembles that of a solid, continued beyond the hrst Bril- 
louin zone. While for n = 128 and a — 0.58 the system is 
still in the gas phase, our PIGS results presented below 
indeed predict a stripe phase at even higher density. 

The dotted lines in Fig. [2] depict the damping limit 
E c (k) above which decay into two excitations of lower 
energy is kinematically allowed, hence excitations be- 
low E c (k) have infinite lifetime corresponding to peaks 
in S(k, E) with zero linewidth. The kinematics of an 
anisotropic dispersion is different from the isotropic case, 
as evidenced e.g. by the lack of a constant Pitaevskii 
plateau. The decay into two rotons is very efficient in 
an isotropic system because of the high density of states 
at the roton energy. For the anisotropic phonon-roton 
dispersion, the roton energy depends on the direction of 
k, thus the roton energies are spread out leading to a 
smoother density of states than in the isotropic limit. 
For example, decay of the maxon in the y-direction is 
not allowed, although its energy is higher than twice the 
roton energy. 

Both the results for S(k, E) and the qualitative sta- 
bility analysis (Fig. [T]) suggest that, as a increases, the 
system develops a preference for long range order in the 
y-direction, until the gas phase becomes unstable at a 
density-dependent critical angle ctQ. The PIGS method 
used to evaluate S(k) in the gas phase can also be used 
to analyze the static properties of a system with long- 
range order as the present one when the homogeneous 
gas is not stable anymore. As in the gas phase, we evalu- 
ate 5(k) because long-range order can be studied by the 
emergence of Bragg peaks. This is indeed what happens 
when a is increased beyond ag. 

Figure [3] summarizes the main PIGS results. The up- 
per left panel shows with black stars the structure factor 
S(k) (shifted up for better visibility) for the isotropic 
(a = 0) system at n = 128, while results along the x- 
and y-directions for n = 64 and a — 0.58 are depicted 
with blue circles and red squares, respectively. S(k, 0) 
and S^O, k) are markedly different in the anisotropic case, 
which is a direct consequence of the anisotropy of the in- 
teraction induced by the non- vanishing tilting angle. Like 
in the isotropic case, the system is in the gas phase ac- 
cording to our stability analysis. We visualize this in the 
upper right panel by a snapshot of one quarter of the sim- 
ulation box corresponding to the (n = 64, a = 0.58) case, 
where each worldline is a different particle. As expected 
for a gas, there is no apparent ordering. 

Results for (n = 128, a = 0.58) and (n = 256, a = 
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FIG. 3. Static structure factor (left panels) and configura- 
tion snapshots (right panels) at different densities and tilt- 
ing angles. S(k, 0) and S(0,k) for n — 64 and a — 0.58 
is shown in the top (blue circles and red squares, respec- 
tively). The black stars show the isotropic S(k) at n — 128 
and a — 0. The middle and bottom panels show the S(k, 
and 5(0, k) and configuration snapshots for n = 128, a = 0.58 
and n = 256, a = 0.61, respectively. 

0.61) are shown in the middle and lower panel, re- 
spectively. As can be seen, the system becomes more 
anisotropic for larger density, the peak in 5(0, k) is more 
pronounced, while the peak in S(k, 0) is less affected. For 
(n = 128, a = 0.58) 5(0, k) is still a smooth function, 
with a peak height that is independent of the number of 
particles N in the simulation. Hence, it is not a Bragg 
peak and the system is still in the gas phase. However, for 
(n = 256 and a — 0.61), the peak in 5(0, k) is orders of 
magnitude larger than the peak in S(k, 0) (notice the log- 
arithmic scale of the bottom panel). The corresponding 
snapshot shows clearly the formation of a stripe phase, 
which according to 5(k) is like a gas in the x-direction 
where the interaction is weak, and a solid along the y-axis 
where the interaction is strong. The peak in 5(0, k) grows 
almost linearly with N, which further supports its inter- 
pretation as a Bragg peak. A second peak of less but still 
significant intensity develops at twice the wave number 
of the first peak. We note that for a = a stripe phase 
has not been observed and the isotropic system remains 
in the gas phase until it solidifies at high density [2] . The 
same conclusion for a fully isotropic interaction has been 
reported recently for a dipolar Fermi gas in 2D [15] . 

We close the discussion showing a 3D-map of the pair 
distribution function g(r) for |r| < L-/2 (£_ is the 
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FIG. 4. Pair distribution function g(x, y) for density n = 128 
and a = 0.58 (left), and for n = 256 and a = 0.61 (right). 

smaller side of the simulation box) in Fig. [4] for the two 
cases (n = 128, a = 0.58) (left panel) and (n = 256, a = 
0.61) (right panel). The stripe phase becomes clearly vis- 
ible in the second case as a plane wave in the y-direction, 
filling the whole simulation box. For (n = 128, a = 0.58) 
oscillations in the y-direction are present only for small 
x and are damped with increasing y, hence g(r) becomes 
isotropic for large |r| and equal to unity, consistent with 
the behavior of a gas. 

Summarizing, we have analyzed the behavior of an 
anisotropic dipolar Bose gas in 2D, using several meth- 
ods: a qualitative stability analysis of the ground state 
based on the HNC-EL method, the exact calculation 
of structural quantities from PIGS Monte Carlo simula- 
tions, and the dynamic structure function in the pair fluc- 
tuation approximation of the dynamic many-body the- 
ory. All results show that for large tilting angle a and 
large density n, there is a quantum phase transition to 
a stripe phase, characterized by long-range order in the 
direction of stronger interaction. The phonon-roton dis- 
persion is very anisotropic in the gas phase close to the 
stripe phase transition, with an almost vanishing roton 
energy in the y-direction. 
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